Order of operations in this file

Load in dplyr, lubridate, and ggplot libraries

library(lubridate)  ## for date/time handling
library(ggplot2)  ## for plotting
library(dplyr)  ## for groupby/summarize analysis

Set a few basic properties for a ggplot theme

cust_theme <- theme(legend.text = element_text(size=12),
                      legend.title = element_text(size=18),
                      axis.text = element_text(size=12),
                      axis.title = element_text(size=18))

Read in the data, handle NA strings

## to run this code, please make sure the three .csv files are in your working directory - thanks!
pressure_temperature <- read.csv("baseline-pressure-temp.csv", skip = 3, colClasses = c("character", "numeric", "numeric", "numeric"), na.strings = c("#N/A", "<NoData>"))

load_power <- read.csv("baseline-load-power.csv",
                       colClasses = c("character", "numeric", "numeric", "numeric"))

exp_period <- read.csv("efficiency.csv", skip = 5,
                       colClasses = c("character", "numeric", "numeric", "numeric", "numeric"),
                       na.strings = c("#N/A", "<NoData>"))

Add a few columns to the dataframes:

Also below, the control series sheets from the original Excel file are merged:

## add a usable date column to each dataset
pressure_temperature$date <- mdy_hm(pressure_temperature$Date...Time)
load_power$date <- mdy_hm(load_power$Date...Time)
exp_period$date <- mdy_hm(exp_period$Date...Time)

## merge the pressure/temp dataset with the load power dataset
control_series <- merge(pressure_temperature, load_power, by="date")

## adding columns to control series:
## add a month column for easy grouping by month and a weekday column for grouping by day of the week
control_series$month <- month(control_series$date)
control_series$weekday <- wday(control_series$date)
exp_period$weekday <- wday(exp_period$date)

## add a column with load in kW instead of tons so the coefficient of performance will be unitless
control_series$load.kW <- control_series$Load..tons. * 3.51685
exp_period$load.kW <- exp_period$Load..tons. * 3.51685

## add a coefficient of performance
control_series$COP <- control_series$load.kW / control_series$Refr.compressor.power..kW.
exp_period$COP <- exp_period$load.kW / exp_period$Refr.compressor.power..kW.

Exploring the Control Year

High wet bulb temperatures can hurt efficiency of the condenser fans. First, let’s see if this was an issue during the January - April period of 2012 (our experimental period spans January through April of 2013).

## How big a problem might high wet bulb temps present?
## identify occurrences when condenser fans are not lowering the discharge pressure to within 1 SD of target
## find standard deviation of discharge pressures over 2012
P_one_sd <- sd(control_series$Discharge.pressure..psig., na.rm=TRUE)
P_mean <- mean(control_series$Discharge.pressure..psig., na.rm=TRUE)

## flag points above 1 SD of discharge pressure
control_series$highP <- cut(control_series$Discharge.pressure..psig., 
                            breaks = c(0, P_mean + P_one_sd, 500), 
                            labels = c("normal", "high"))

## add a month name column to label the facet plot
control_series$month_name <- factor(month.name[control_series$month], 
                                    levels= c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"))


## how many "high pressure" periods by month and what fraction do they represent:
P_by_month <- group_by(control_series, month, highP) %>%
                summarize(n = n()) %>%
                mutate(fraction = n / sum(n))
high_P_by_month <- filter(P_by_month, highP == "high")

## plot discharge pressure vs. wet bulb temp, highlighting those points more than 1 SD above target
g <- ggplot(control_series, aes(x=Twb..F., y=Discharge.pressure..psig., color=highP)) + 
    geom_point(alpha=.12,na.rm = TRUE) +
    scale_color_manual("Pressure Status\n", labels=c("Normal", "High"), values=c("#000080", "red"))
    coord_cartesian(xlim = c(-10, 80), ylim = c(70,180))
## <ggproto object: Class CoordCartesian, Coord>
##     aspect: function
##     distance: function
##     expand: TRUE
##     is_linear: function
##     labels: function
##     limits: list
##     range: function
##     render_axis_h: function
##     render_axis_v: function
##     render_bg: function
##     render_fg: function
##     train: function
##     transform: function
##     super:  <ggproto object: Class CoordCartesian, Coord>
g + facet_wrap( ~ month_name, ncol=3) + 
    cust_theme + 
    xlab(bquote('Wet Bulb Temperature ('^o*F*')')) + 
    ylab("Actual Discharge Pressure (psi)") +
    guides(colour = guide_legend(override.aes = list(alpha = 1)))

Points highlighted on the above graph are those that are more than one standard deviation above the target (145 psi).

Here’s the fraction of each of these months that see pressures that are too high (more than 1 SD above target):

## plot fraction of "high" pressure readings by month
ggplot(high_P_by_month, aes(x=factor(month), y=fraction)) + 
    geom_bar(stat="identity") +
    geom_text(aes(color="white", label=round(fraction, digits=2)), vjust=1.5) +
    cust_theme +
    theme(legend.position="none") +
    xlab("Month") +
    ylab("Fraction of All Readings That Are High Pressure")

So, while wet bulb temperatures may cause compressor pressure problems in the summer, we should have minimal problems looking at an experiment that starts in January and ends in mid-April.

Suffering COP

Even though our experimental period should be relatively unaffected by high wet bulb temperatures, it would be interesting to know how much the compressor’s efficiency suffers with high wet bulb temperatures.

Below, I’ve plotted the compressor’s coefficient of performance by month and “load category”. A couple trends worth noting: * The compressor seems to get more efficient with greater load. * Compressor COP standard deviation drops with greater load. * With high load and high wet bulb temperature, the COP also drops.

## how does the COP vary with wet bulb temperature at different loads?
## cut the load index into different bins
control_series$load_category <- cut(control_series$load.kW,
                                    breaks = c(-1, 1, 200, 800, 1200, 2200),
                                    labels = c("No Load", "Very Low", "Low", "Medium", "High"))

## plot of all months and load categories (very low < 200kW < low < 800kW < medium < 1200kW < high < 2200kW)
h <- ggplot(filter(control_series, load_category != "No Load"), aes(x=Twb..F., y=COP)) +
    geom_point(aes(alpha = .12, color = Twb..F.), na.rm = TRUE) + 
    scale_color_gradient(low="blue", high="red")
h + facet_grid(month_name ~ load_category) + geom_rug(position="jitter", size = 0.1) +
    cust_theme +
    xlab(bquote('Wet Bulb Temperature ('^o*F*')')) +
    ylab("Coefficient of Performance") +
    ggtitle("Load Category") + 
    theme(legend.position="none")

Here’s an example of a summer month that sees a real drop in COP at high load:

## example plot of COP suffering because condenser fans aren't removing enough heat before refrigerant returns to the compressor
## July 2012
ggplot(filter(control_series, load_category == "High", month == 7), aes(x=Twb..F., y=COP)) +
    geom_point(aes(alpha = .12, color = Twb..F.), na.rm = TRUE) + 
    scale_color_gradient(limits=c(50, 77),low="blue", high="red") +
    coord_cartesian(xlim = c(48, 80)) +
    cust_theme + 
    theme(legend.position="none") +
    ggtitle("July 2012") +
    xlab(bquote('Wet Bulb Temperature ('^o*F*')')) +
    ylab("Coefficient of Performance")

Because the experimental periods are short (two of them are only about 10 days), I thought it might also be worthwhile to make sure I didn’t bias those experiments by weekday. If one 10 day period has two full weekends whereas another has only one weekend, that might make a big difference. Let’s check to see if there’s a significant difference in load profile by weekday:

## how much does the load depend on the day of the week?
load_by_month_wday <- group_by(control_series, month_name, weekday, load_category) %>%
    summarize(n = n()) %>%
    mutate(wday_fraction = n / sum(n))

i <- ggplot(load_by_month_wday, aes(x=weekday, y=wday_fraction, fill=load_category)) + geom_bar(stat="identity") + scale_fill_hue(c=45, l=80)
i + facet_wrap(~ month_name) +
    cust_theme +
    ylab("Load Fraction") +
    xlab("Day of Week") +
    labs(fill="Load") +
    scale_x_discrete(name="Day of Week", limits=c("S", "M","T","W","R","F","Sa"))

It’s clear that the load profile changes dramatically for day of the week. This refrigerator clearly doesn’t see the same load on the weekends as it does during the week. I’ll have to compensate for this during the analysis of the experimental period.

Conclusions from exploring control data:

  • The experimental period from January through April 2013 should see minimal interference from high wet bulb temps (though it will be hard to predict the different pressures’ behaviors with high wet bulb temperatures).
  • We’ll have to find representative coefficients of performance for each experimental pressure and load category, then extrapolate those to year-long performance assuming a similar load profile to 2012.

The Experiment

First, define when each period of the experiment started and finished

On days when the target pressure was changed, it’s not clear exactly when the change was made, so I’ve thrown those dates out entirely from the analysis in this section.

## define intervals for the experimental data
exp_period$p_label <- cut(exp_period$date, breaks= as.POSIXct(c("2012-12-31", "2013-02-03", "2013-02-05", "2013-02-19", 
                                                                "2013-02-21", "2013-03-06", "2013-03-08", "2013-05-21")),
                                            labels = c("p145", "t1", "p140", "t2", "p135", "t3", "p130")) 

exp_period$load_category <- cut(exp_period$load.kW,
                                    breaks = c(-1, 1, 200, 800, 1200, 2200),
                                    labels = c("No Load", "Very Low", "Low", "Medium", "High"))

Check load similarity among experimental periods

Recall from the earlier section that refrigeration load varies dramatically during the week and the test periods are relatively short. If we are going to compare the COPs from each experiment directly, the load profile should be similar among all experiments…

## lets make sure the load profile is similar for all the periods in question
load_by_period <- group_by(exp_period, p_label, load_category) %>%
    summarize(n = n()) %>%
    mutate(load_fraction = n / sum(n)) %>%
    filter(p_label != "t1" & p_label != "t2" & p_label != "t3")

## the load profiles reveal that the comparison isn't really all that fair:
ggplot(load_by_period, aes(x=p_label, y=load_fraction, fill=load_category)) + 
    geom_bar(stat="identity") + 
    scale_fill_hue(c=45, l=80) +
    cust_theme +
    ylab("Load Fraction") +
    xlab("Experiment Period (pressure setting (psi))") +
    labs(fill="Load")

… but it (the load profile) varies considerably from one experiment to the next. The approach from here will be to compute an average COP for each load category (Very Low, Low, Medium, and High) and target compressor pressure (130, 135, 140, and 145 psi) and use those COPs combined with 2012’s year-long load category profile to create an overall net COP for each pressure setting.

First, let’s calculate the COP for each compressor pressure and load category, weighting the average by how much energy the compressor is consuming. A COP taken while the compressor is consuming 50 units of energy is weighted 50 times more heavily than a COP calculated when the refrigerator is consuming 1 unit of energy:

## let's compare COPs for each load group independently among the 4 experimental groups
## COP averages weighted by load
weighted_means <- group_by(exp_period, p_label, load_category) %>%
    summarize(n = n(), sum_load = sum(load.kW), sum_compressor_power = sum(Refr.compressor.power..kW.)) %>%
    mutate(w_mean = sum_load/sum_compressor_power) %>%
    filter(p_label != "t1" & p_label != "t2" & p_label != "t3" & load_category != "No Load")

ggplot(weighted_means, aes(x=load_category, y=w_mean, fill=p_label)) +
    geom_bar(stat="identity",position="dodge") + 
    scale_fill_hue(c=45, l=80) +
    cust_theme +
    ylab("Compressor COP") +
    xlab("Load Category") +
    labs(fill="Target Pressure")

While there’s no real trend in COP at the “Low” and “Very Low” load categories, there’s a very clear trend at “Medium” and “High” loads. Great – the experiment was at least a partial success!

One more step before we get to the bottom line: we need to find typical overall loads for an entire year – we’ll use the 2012 data for that:

## now, assuming a load distribution like 2012's, let's estimate an average COP for each target pressure setting
## first, what's the load distribution like for 2012?
year_long_load <- sum(control_series$load.kW)
control_load_fractions <- group_by(control_series, load_category) %>%
    summarize(number = n(), load_fraction = sum(load.kW)/year_long_load)

## now let's assume those load fractions will be similar for 2013 and calculate expected COPs for each load category - pressure combination:
year_long_predicted_COPs <- merge(weighted_means, control_load_fractions, by="load_category") %>%
    mutate(COP_contribution = load_fraction * w_mean) %>%
    group_by(p_label) %>%
    summarize(expected_COP = sum(COP_contribution))
print(year_long_predicted_COPs)
## Source: local data frame [4 x 2]
## 
##   p_label expected_COP
##    (fctr)        (dbl)
## 1    p145     4.907121
## 2    p140     5.202618
## 3    p135     5.288401
## 4    p130     5.363794
## plot the expected year long COPs
ggplot(year_long_predicted_COPs, aes(x=p_label, y=expected_COP, fill=p_label)) +
    geom_bar(stat="identity", position="dodge") +
    scale_fill_hue(c=45, l=80) +
    coord_cartesian(ylim = c(4.5,5.5)) +
    cust_theme +
    ylab("Predicted Annual Net Compressor COP") +
    xlab("Compressor Target Pressure") +
    theme(legend.position="none")

Great! There are the predicted, load-weighted, experiment-produced, year-long COPs for each compressor target pressure. And good news: it looks like the experiment has produced exactly the results we’d hoped for – the lowest compressor pressure should be the most efficient.

Let’s see if we can turn that into a year-long savings (assuming $0.15 per kWh). # The Important Part:

## energy projected for 365 days of 2012 -- control series is missing a few dates
correction_factor <- (365*24*4)/nrow(control_series)

total_year_heat_pumped <- sum(control_series$load.kW) * .25 * correction_factor # 15 minute data requires .25 correction factor
projected_savings <- mutate(year_long_predicted_COPs, annual_predicted_energy = total_year_heat_pumped / expected_COP,  Ecost_per_year = annual_predicted_energy * .15) # assuming $0.15 per kWh

## simple bar graph of the savings
ggplot(projected_savings, aes(x=p_label, y=Ecost_per_year, fill=p_label)) + 
    geom_bar(stat="identity") +
    coord_cartesian(ylim=c(200000, 250000)) +
    scale_fill_hue(c=45, l=80) +
    cust_theme +
    ylab("Predicted Annual Cost to Run Compressors (USD)") +
    xlab("Compressor Target Pressure") +
    theme(legend.position="none")

Woohoo! Switching from a compressor pressure of 145 psi to 130 psi finds a savings of about $21k or ~9% of the total energy used by the compressor over the course of the year!

The little caveat

The above savings estimate is almost certainly a bit high. The condenser fans are going to have to work harder at the lower compressor pressures to lower the temperature of the refrigerant, offsetting part of the savings. Also at lower compressor pressures, the condenser fans will have trouble doing their jobs at a lower wet-bulb temperature than at higher compressor pressures (where we were already seeing some problems with the compressor!).

We were given the load of the entire rest of the facility to try to detect if there was any increased use in condenser fan usage, but the data appears to be too noisy to detect this increased condenser fan load. Perhaps if the summer were repeated in the summer the condenser fan load would be more apparent. Even better, get the exact power consumed by the condenser fans so a net condenser fan plus compressor COP could be calculated.

Here’s a quick attempt I made to see if I could detect an increase in condenser fan usage within the “other” energy category. If the condenser fans are running more at lower compressor target pressures and all else is equal, the “other” power should be highest for the lowest compressor target pressure, especially under high load conditions. If that is happening, it’s not very clear from the data:

## first, get a "non-compressor" power consumption column
exp_period$other_power <- (exp_period$Utility.power..kW. - exp_period$Refr.compressor.power..kW.)

## now, we should see an increase in average "non-compressor" power at lower compressor target pressures
non_compressor_power_means <- group_by(exp_period, p_label, load_category) %>%
    summarize(cnpm = mean(other_power)) %>%
    filter(p_label != "t1" & p_label != "t2" & p_label != "t3" & load_category != "No Load")

ggplot(non_compressor_power_means, aes(x=load_category, y=cnpm, fill=p_label)) +
    geom_bar(stat="identity", position = "dodge") +
    scale_fill_hue(c=45, l=80) +
    cust_theme +
    ylab("Non-compressor Power Means (kW)") +
    xlab("Compressor Load Category") +
    labs(fill="Target Pressure")

Thanks for reading and the candid advice! - Mike